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1.0  INTRODUCTION 


Muzzle  blast  and  flash  are  inherent  to  ejection  of  a  chemically-propelled 
projectile  from  a  gun  barrel.  Blast  is  the  result  of  release  of  high-pressure 
gases  and  the  ensuing  over-pressure  of  the  surrounding  atmosphere.  Flash  is 
the  radiation  from  the  expelled  propellant  gases,  which  is  frequently  augmented 
by  afterburning  of  the  hydrogen-rich  propellent  gases  in  the  atmosphere.  Blast 
and  flash,  as  well  as  the  dispersal  of  propellant  gases  and  combustion 
products,  represent  environmental  hazards  to  the  gunner  and  other  personnel  or 
machinery  in  the  area  of  the  gun  platform.  Additionally,  they  are  functionally 
undesirable  in  that  they  reveal  the  gun  location.  Consequently,  considerable 
effort  has  been  devoted  to  development  of  mechanical  devices  which  are  intended 
to  suppress  blast  and  flash  or  to  deflect  muzzle  gases  away  from  sensitive 
areas  (e.g.,  aircraft  engine  inlets).  Although  blast,  flash  and  the 
distribution  of  propellent  exhaust  gases  are  closely  interrelated  phenomena, 
muzzle  devices  which  are  successful  in  achieving  one  objective  (e.g.,  exhaust 
gas  diverters  and  recoil  attenuation  devices)  are  often  counterproductive  in 
other  areas  (i.e.,  flash  suppression). 


The  coupling  between  muzzle  blast  and  flash  is  evident  in  the 
amplification  of  blast  overpressure  by  propellant  afterburning,  in  the 
augmentation  of  flash  by  muzzle  brakes  (which  divert  the  exhaust  gases),  and  in 
flash  suppression  by  conical  exhaust  nozzles  (which  focus  the  exhaust  gases). 
Blast  and  flash  therefore  clearly  involve  an  interplay  between  mechanical 


design  of  the  gun  and  chemical  composition  of  the  propellant.  A  more  thorough 


discussion  of  the 
studies  by  Keller 


phenomenology  and 
*  ^  and  Schmidt^ 


recent  experience  can  be  found  in  the 
as  well  as  in  the  substantial  body  of  prior 


work  cited  therein. 


The  present  work  concentrates  on  development  of  a  model  of  muzzle  flash 
which  embodies  the  essential  thermo -chemistry  and  fluid-dynamics  of  the  process 
while  retaining  computational  simplicity  appropriate  to  exploratory  studies  and 
engineering  analyses.  The  essential  features  of  the  phenomenon  are  considered 
to  be  unsteady  mixing  of  the  propellant  exhaust  gases  with  shock-heated  air  and 
concurrent  chemical  reaction  in  a  flowfield  dominated  by  an  expanding  blast 
wave.  Emphasis  has  been  placed  on  identification  and  representation  of  the 
primary  physical  processes  rather  than  detailed  numerical  analysis  or 
comprehensive  comparative  studies. 


1  Keller,  G.E. ,  "Secondary  Muzzle  Flash  and  Blast  of  the  British  81-mm, 
L16A2,  Mortar"  ARBRL-MR-03117 ,  July  1981. 

^  Keller,  G.E.,  "The  Effect  of  Propellant  Composition  on  Secondary 
Muzzle  Blast  Overpressure"  ARBRL-MR-03270,  April  1983. 

^  Schmidt,  E.M.,  "Secondary  Combustion  in  Gun  Exhaust  Flows"  ARBRL-TR- 
02373,  October  1981.  — 
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The  basic  fluid-dynamical  model  draws  on  the  muzzle  blast  analysis  of 
Erdos  and  DelGuidice^,  which  describes  the  unsteady  flowfield  along  the  gun 
axis  from  the  muzzle  to  the  blast  wave.  This  model  assumes,  however,  that  the 
propellant  gases  are  effectively  inert  after  expulsion  from  the  muzzle  and 


neglects  any  mixing  along  their  interface  with  the  atmosphere.  Both  these 
assumptions  are  removed  in  the  present  model.  The  chemical-kinetic  mechanisms 


governing  the  reaction  of  the  propellant  gases  with  air,  including  the  effects 
of  potassium  additives,  are  drawn  from  the  work  of  Yousef ian^,^.  This  work 


pays  careful  attention  to  the  complex  kinetics  of  the  afterburning  reactions 


and  the  potential  for  flash  suppression  by  chemical  additives.  However,  the 


fluid-dynamics  are  modeled  by  steady-state  mixing  with  constant  boundary 
conditions,  which  unfortunately  ignores  essential  transients  in  pressure  and 
variations  in  boundary  conditions  on  the  mixing  layer  due  to  the  blast  wave.* 


The  present  work  also  draws  on  the  study  of  secondary  combustion  by 
Schmidt-*.  Schmidt  recognizes  the  essentially  unsteady  nature  of  the 
phenomenon,  but  attempts  to  apply  the  unsteady  conditions  to  empirical  flash 
criteria  without  modeling  the  rate  of  mixing  or  of  reaction.  His  study  is  also 
important  in  that  it  identifies  the  mechanism  by  which  muzzle  brakes  tend  to 
augment  flash,  namely  by  the  generation  of  multiple  shocks.  The  present  work 
does  not  attempt  to  explicitly  model  muzzle  devices,  such  as  brakes.  However 
it  does  allow  the  effects  of  such  devices  to  be  modeled  through  the  fluid- 
dynamic  boundary  conditions  placed  on  the  mixing  layer,  in  the  same  fashion  as 
Schmidt  considers  their  effect  on  the  empirical  flash  criteria. 


Finally,  the  evaluation  of  available  muzzle  flash  codes  performed  by 
Keller^  is  pointed  out  since  it  reinforces  the  requirement  for  a  model  such  as 
proposed  herein. 


A  Erdos,  J.  and  DelGuidice  P.,  "Calculation  of  Muzzle  Blast  Flowfields" 
AIAA  Journal,  Vol.  13,  No.  8,  August  1975  pp.  1048-1056. 

^  Yousefian,  V.,  "Muzzle  Flash  Onset"  ARBRL-CR-00477 ,  Feb.  1982. 

®  Yousefian,  V.,  "Muzzle  Flash  Onset:  An  Algebraic  Criterion  and 
Further  Validation  of  the  Muzzle  Exhaust  Flow  Field  Model"  ARBRL-CR-00506, 
March  1983. 

^  Keller,  G.E.,  "An  Evaluation  of  Muzzle  Flash  Prediction  Codes"  ARBRL- 
MR-03318,  Nov.  1983. 
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DISCUSSION  OF  THE  MODEL 


The  essential  features  of  a  muzzle  exhaust  flowfield  are  sketched  in 
Figure  (1).  The  propellant  gases  expand  supersonically  upon  release  from  the 
gun  barrel.  In  a  vacuum,  the  expansion  process  would  approach  a  spherical 
source  flow  with  it's  origin  at  the  muzzle.  The  propellant  gases  would  form, 
in  effect,  an  expanding  spherical  "balloon".  In  the  atmosphere,  the  spherical 
expansion  is  contained  by  a  finite  back-pressure  exerted  by  the  surrounding 
air.  This  boundary  condition  leads  to  the  barrel  shocks  which  terminate  the 
near-spherical  expansion  and  turn  the  propellant  gases  back  to  a  near-axial 
flow  direction.  The  expanding  front  of  propellant  gases  still  forms  a 
"balloon"  of  sorts  which  drives  a  strong  shock  wave,  i.e.,  a  "blast  wave",  into 
the  atmosphere  ahead  of  it.  The  pressure  field  behind  the  blast  wave  causes 
the  supersonic  propellant  gases  to  compress  to  a  matching  pressure  across 
another  strong  shock,  termed  the  "Mach  disc".  The  subsonic  flowfield  between 
the  Mach  disc  and  the  blast  wave  again  resembles  a  spherical  expansion,  a  fact 
that  allowed  Erdos  and  DelGuidice^  to  model  this  region  as  a  one-dimensional 
unsteady  flow  with  spherical  symmetry. 

The  supersonic  portion  of  the  propellant  expansion  is  quasi-steady;  i.e., 
it  adjusts  almost  instantly  to  the  time-varying  conditions  at  the  muzzle  as  the 
propellant  gases  empty  from  the  gun.  The  position  of  the  Mach  disc  varies 
continuously,  but  upstream  of  the  disc  the  flowfield  resembles  a  steady,  under¬ 
expanded  rocket  plume.  Mixing  of  the  propellant  gases  with  air  will  occur 
along  the  lateral  boundaries  of  the  plume,  and  since  the  propellant  gases  are 
typically  rich  in  hydrogen  and  carbon  monoxide,  some  afterburning  in  this 
mixing  layer  may  be  expected.  However,  since  the  air  side  of  the  layer  is 
essentially  quiescent  while  the  propellant  gases  are  moving  at  high  supersonic 
speeds,  and  since  the  total  (stagnation)  temperature  of  the  air  is  much  lower 
than  that  of  the  propellant  gases,  the  maximum  static  temperature  in  the  mixing 
layer  can  be  estimated  by  analogy  with  a  hypersonic  cold-wall  boundary  layer®. 
Accordingly,  the  maximum  static  temperature  will  be  about  one-quarter  the  total 
temperature  of  the  propellant  gases,  and  this  temperature  will  occur  where  the 
velocity  is  about  one-half  the  propellant  velocity  and  the  gas  mixture  will 
consist  of  about  one-half  air  and  one-half  propellant  gases.  The  mij  are  ratio 
and  temperature  would  therefore  favor  ignition,  but  the  velocity  may  inhibit 
ignition  by  allowing  less  time  at  this  condition  than  the  required  ignition 
delay  time.  Characteristic  gas  residence  time  in  this  layer  will  be  a  small 
fraction  of  the  gun  emptying  time. 

The  subsonic  portion  of  the  propellant  gas  expansion  i.e.,  the  portion 
between  the  Mach  disc  and  the  propellant-air  interface,  is  decidedly  non¬ 
steady.  Since  the  Mach  number  of  the  flow  is  low,  the  static  temperature  of 
the  propellant  gases  approaches  the  total  temperature.  Furthermore,  the 
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velocity  is  low  so  the  residence  time  in  this  layer  will  approach  the  gun 
emptying  time.  Consequently,  afterburning  is  far  more  likely  in  the  mixing 
layer  that  develops  at  the  propellant-air  interface  in  this  region  (i.e.,  near 
the  gun  axis)  than  along  the  boundaries  of  the  supersonic  plume. 

The  static  temperature  of  the  propellant  gas  and  of  the  air  at  their 
interface  varies  in  time  as  the  pressure  drops.  However,  the  entropy  of  the 
propellant  gas  and  of  the  air  at  any  point  in  the  flowfield  defines  the 
thermodynamic  relationship  between  the  temperature  and  pressure.  The  time 
variation  of  entropy  is  governed  by  the  energy  equation  which  in  the  absence  of 
chemical  reactions,  heat  conduction,  diffusion,  etc.,  simply  states  that  the 
entropy  is  convected  along  streamlines  (particle  paths)  except  where  the 
streamlines  cross  shock  waves.  The  entropies  of  the  two  sides  of  the  interface 
are  therefore  determined  by  the  initial  conditions  at  release  of  the  propellant 
gas.  The  blast  wave  begins  as  a  strong  shock  wave  and  maintains  its  strength 
as  the  cumulative  amount  of  energy  released  by  the  escaping  propellant  gases 
increases.  The  strength  of  the  blast  wave  begins  to  decay,  however,  as  it 
acquires  a  spherical  shape  and  grows  in  radius.  The  cumulative  energy  released 
by  the  propellant  gases  also  begins  to  asymptote  to  a  constant  as  the  gun 
barrel  empties,  and  the  blast  wave  decay  rate  transitions  from  the  power  law 
associated  with  a  constant  rate  of  energy  deposition  to  a  constant  total  energy 
release  (i.e.,  a  point  explosion).  The  net  result  is  that  the  entropy  of  the 
shock-heated  air  is  initially  about  a  constant  near  the  propellant-air 
interface  and  then  steadily  decreases  with  increasing  distance  from  the 
interface. 

The  entropy  of  the  propellant  gas,  on  the  other  hand,  consists  of  the 
entropy  generated  by  the  propellant  combustion  process  in  the  breech-end  of  the 
gun  barrel  plus  the  increment  added  as  the  gas  passes  across  the  plume  Mach 
disc.  The  entropy  layer  on  the  propellant  side  of  the  interface  therefore  also 
exhibits  a  spatial  variation  associated  with  the  initial  increase  in  strength 
of  the  Mach  disc,  but  as  the  Mach  disc  later  decays  in  strength  and  collapses 
the  propellant-side  entropy  asymptotes  to  the  combustion  generated  value.  The 
propellant-side  flowfield  therefore  exhibits  a  significant  increase  in  entropy 
with  increasing  distance  from  the  interface,  a  maximum  entropy  at  some  point, 
and  a  subsequent  decrease  and  asymptotic  approach  to  the  entropy  of  the  gas 
venting  from  the  muzzle. 

The  significance  of  the  propellant-side  and  air-side  entropy  layers  is 
that  the  shock-processed  air  is  hot  as  it  accumulates  near  the  interface,  and 
remains  relatively  warm  as  the  pressure  drops  and  volume  expands.  However,  the 
hot  propellant  gas  leaving  the  muzzle  is  further  heated  by  the  Mach  disc,  and 
as  the  pressure  drops  a  temperature  distribution  resembling  the  entropy  layer 
develops,  namely:  warm  at  the  interface,  hotter  proceeding  away  from  the 
interface  to  a  maximum  temperature,  and  then  cooling  again  and  asymptoting  to 
the  temperature  of  the  gas  venting  through  an  isentropic  subsonic  expansion  to 
the  atmosphere.  Consequently,  as  the  gases  mix  along  the  interface,  the  mixing 
layer  will  initially  entrain  relatively  hot  propellant  gas  and  air.  (The  term 
"relative"  is  used  to  denote  the  dependence  on  pressure. )  As  more  propellant 
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gas  is  entrained,  its  relative  temperature  will  increase,  accelerating  the  rate 
at  which  afterburning  reactions  between  the  air  and  the  unbumed  propellant  gas 
proceed.  The  relative  temperature  of  the  entrained  propellant  gas  will  reach  a 
maximum  and  then  slowly  decay  toward  an  asymptotic  final  temperature. 


Reactions  between  the  unburned  propellant  gas  and  air  are  initiated  in  the 
mixing  layer  at  very  high  pressure,  but  as  the  pressure  drops  the  molecular 
collision  frequency  drops  in  direct  proportion  and  the  speed  of  the  chain¬ 
forming  ignition  reactions  may  not  be  sufficient  to  sustain  combustion. 
However,  the  ignition- sustaining  reactions  will  be  promoted  by  the  entrainment 
of  increasingly  hotter  propellant  gas,  which  will  offset  the  effect  of  dropping 
pressure  to  some  extent,  since  the  reaction  rates  are  typically  exponentially 
proportional  to  temperature.  If  the  ignition  reactions  become  self-sustaining 
during  the  period  of  entrainment  of  shock  heated  propellant  gas,  flash  occurs. 
If  they  have  not  become  self-sustaining  by  the  time  the  entrained  propellant 
gas  has  reached  it's  asymptotic  decay  period,  there  will  be  no  flash  due  to 
afterburning. 


Clearly  then,  the  occurrence  of  flash  is  critically  dependent  on  the 
competing  mechanisms  of  falling  pressure  and  entrainment  of  variable 
temperature  propellant  gas  which  are  superimposed  on  the  mixing  process.  The 
present  model  has  been  specifically  formulated  to  address  this  phenemonology. 


The  role  of  chemical  additives  in  the  propellant  is  simply  to  interfere 
with  the  ignition  process  by  tying  up  the  free-radicals  needed  to  sustain 
combustion.  The  amount  of  additive  needed  will  depend  on  the  transient  fluid- 
mechanical  characteristics  of  the  specific  muzzle  flowfield.  The  present  model 
is  also  intended  to  provide  a  capability  to  rapidly  assess  the  effectiveness  of 
chemical  additives. 


As  the  preceding  discussion  of  the  phenomology  should  indicate,  the 
transient  characteristics  of  the  flowfield,  particularly  the  pressure  impressed 
on  the  mixing  layer  and  the  variation  in  temperature  (i.e.,  entropy)  of  the 
entrained  gas,  are  considered  essential  to  evaluation  of  the  chemical  processes 
leading  to  flash.  These  processes  occur  on  a  time  scale  comparable  to  the 
characteristic  time  for  the  flowfield  transients  and  must  therefore  be  modeled 
as  "finite  rate"  processes.  Unfortunately,  the  cost  of  carrying  out  finite- 
rate  chemical  calculations  can  escalate  very  quickly.  The  fluid  dynamics 
require  solution  of  three  equations:  conservation  of  mass,  momentum  and 
energy.  Finite-rate  chemistry  adds  one  more  equation  for  each  chemical  species 
considered,  and  requires  evaluation  of  complex  reaction  rate  expressions  for 
every  reaction  in  which  the  considered  species  participates.  Consequently 
there  is  a  strong  fiscal  motivation  for  keeping  the  number  of  spatial  locations 
at  which  the  species  equations  must  be  solved  to  a  minimum.  The  bare  minimum 
is,  of  course,  one  point,  e.g.,  the  interface. 


A  two-layer  model  of  the  mixing  layer  has  been  adopted  in  which  the 
inviscid  interface  divides  the  two  layers.  The  governing  equations  for 
conservation  of  mass  and  energy  are  analytically  integrated  from  the  interface 
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outward  to  the  "edge"  or  "front"  of  each  layer,  leading  to  a  pair  of  ordir.ar. 
differential  equations  for  the  thermal  thickness  of  each  layer.  The  momer*  _j:, 
equation  is  eliminated  from  present  consideration  by  restricting  attention  • 
the  axis  of  symmetry.  The  integral  species  conservation  equation  is 
automatically  satisfied  by  assuming  the  thickness  of  the  species  diffusion 
layer  is  equal  to  the  thermal  thickness. 

Thus,  the  present  model  only  requires  solution  of  the  governing  equations 
at  the  interface,  plus  the  two  equations  for  the  thermal  thickness  of  each 
layer. 

A  detailed  description  of  the  analytical  formulation  of  the  two- layer 
model  is  presented  in  the  following  section. 
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3.0  DESCRIPTION  OF  THE  ANALYTICAL  FORMULATION 


Assumptions  and  Basic  Equations 

As  previously  depicted  in  Figure  (1),  a  mixing  layer  develops  along  the 
interface  between  the  propellant  gas  and  the  air.  Within  this  mixing  layer, 
hot  propellant  gas  (shock-heated  by  the  Mach  disc)  is  brought  into  contact  with 
hot  oxygen  (in  the  air  which  has  been  shock-heated  by  the  blast  wave).  The 
mixing  layer  can  be  viewed  as  two  sublayers:  one  which  extends  outward  from 
the  inviscid  interface  between  the  propellant  gas  and  air  and  one  which  extends 
inward  from  it.  The  inviscid  interface  provides  a  convenient  frame  of 
reference  for  the  mixing  layer. 

The  present  model  is  restricted  to  development  of  a  planar  mixing  layer  at 
the  gun  axis.  In  other  words,  the  radius  of  curvature  of  the  inviscid 
interface  is  assumed  to  be  much  larger  than  the  thickness  of  the  mixing  layer. 
Furthermore,  forces  due  to  acceleration  of  the  interface  are  neglected.  The 
latter  assumption  is  perhaps  suspect  at  early  times,  but  the  interface  location 
becomes  virtually  steady  by  the  time  secondary  combustion  is  usually  observed. 
The  former  assumption  may  also  be  questioned  for  later  times  when  the  mixing 
layer  becomes  very  thick.  Nevertheless,  the  present  model  should  adequately 
model  the  phenomena  of  interest,  and  provide  a  basis  for  further  modeling  as 
needed  to  refine  the  predictive  capabililty. 

The  inviscid  interface  is  therefore  used  as  the  base  of  a  two- 
dimensional,  rectangular  coordinate  system:  x  is  the  distance  along  the 
interface  and  y  is  the  distance  normal  to  it. 

Modeling  of  the  diffusive  processes  poses  a  somewhat  perplexing  dilemma. 
Shadowgraphs  of  muzzle  blast  flowfields  such  as  those  taken  by  Schmidt  and 
Shear9  show  highly  turbulent  structures  in  the  mixing  layers.  However,  the 
propellant-air  interface  near  the  gun  axis  is  not  visible.  No  inviscid 
velocity  jump  exists  across  the  interface  at  the  axis,  and  hence  the  usual 
mechanism  for  generating  turbulent  motion  is  absent  at  that  point.  The 
interface  cannot  support  a  pressure  jump,  and  therefore  the  pressure  gradients 
along  the  interface  must  be  identical.  Consequently,  differences  in  the 
velocity  gradients  must  be  associated  with  the  differences  in  density  across 
the  interface.  Hence,  velocity  jumps  that  develop  away  from  the  axis  due  to 
the  pressure  gradient  are  expected  to  be  relatively  small.  Again,  there  seems 
to  be  an  absence  of  strong  driving  mechanisms  for  turbulence  generation  near 
the  axis.  By  contrast,  the  inviscid  velocity  jumps  across  the  lateral 
boundaries  of  the  plume  are  large,  and  the  observed  turbulence  in  this  area  is 
clearly  consistent  with  a  shear-driven  mechanism. 


Schmidt,  E.  and  Shear,  D.  "Flowfield  About  the  Muzzle  of  an  M16 
Rifle"  BRL  Report  1692,  Jan.,  1974.  Also,  AIAA  Journal,  Vol.  13,  PP.  1088- 
1093,  Aug.,  1975. 
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Acceleration  of  an  interface  between  fluids  of  differing  density  in  the 
direction  of  the  lower  density  fluid  gives  rise  to  the  classical  Rayleigh- 
Taylor  instability  of  the  interface.  The  observed  broadening  of  the  interface 
and  turbulent  mixing  of  driver  and  driven  gases  in  a  shock  tube  has  been 
attributed  to  this  mechanism^  .H.12,13.  The  interface  in  a  shock  tube  is 
proposed  herein  as  an  possible  model  for  the  muzzle  flowfield  interface  at  the 
axis  of  symmetry. 

In  view  of  the  lack  of  a  well-defined  basis  for  modeling  the  diffusive 
mechanism  for  the  mixing  layer  that  develops  along  the  interface,  three  models 
have  been  postulated.  The  first  is  simply  a  laminar  viscosity,  thermal 
conductivity  and  diffusivity  model.  The  second  is  a  heuristic  thermal 
conductivity  model  based  on  a  classical  eddy  viscosity-type  formulation.  The 
third  is  a  two-equation  turbulent  thermal  energy  model  developed  by  analogy 
with  the  K-e  model  for  the  turbulent  kinetic  energy  equation.  In  the  latter 
two  models,  the  effective  Lewis  and  Prandtl  numbers  are  assumed  to  be  unity. 

The  propellant  combustion  products  are  described  as  a  reacting  mixture  of 
chemical  species,  each  behaving  as  a  perfect  gas  (in  vibrational  equilibrium). 
Similarly,  air  is  described  as  a  mixture  of  molecular  and  atomic  nitrogen  and 
oxygen,  and  nitric  oxide  (also  in  vibrational  equilibrium).  Fick's  law,  with  a 
single  diffusion  coefficient  for  all  combinations  of  species,  is  used  tc 
describe  mass  diffusion. 


The  thin-layer  approximations  to  the  Navier-Stokes  equations  will  be 
assumed  to  apply,  i.e.: 


32/3y2  >>  32/3x2  and  p  =  p(x,t). 


The  following  set  of  governing  equations  thereby  obtains  (expressed  here  in 
conservative  form) : 


Richtmeyer,  R.  D. ,  "Taylor  Instability  in  Shock  Acceleration  of 
Compressible  Fluids,"  Communications  on  Pure  and  Applied  Mathematics,  Vol. 
13,1960,  pp.  297-319. 

H  Levin,  M.  A.,  "Turbulent  Mixing  at  the  Contact  Surface  in  a  Driven 
Shock-Wave,"  The  Physics  of  Fluids,  Vol.  13,  No.  5,  1970,  pp.  1166-1171. 

12  Andronov,  V.  A.  et  al.,  "Turbulent  Mixing  at  Contact  Surface 
Accelerated  by  Shock-Waves,"  Soviet  Physics- JETP,  Vol.  44,  1976,  pp.  424-427. 

13  Houas,  L.,  Brun,  R.  and  Hanana,  M. ,  "Experimental  Investigation  of 
Shock-Interface  Interactions,"  AIAA  J,  Vol.  24,  No.  8,  August  1986,  pp.  1254- 
1255. 
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3t  3x  3y  3y  kcp  3y;  3t 


3pa^  3pua^  3pva^ 
3t  +  3x  +  3y 


a  3olj 

fe  (pnr  jji)  +  4i 


with  the  following  equations  of  state: 

p  »  pRQTZCaj/MWj) 

h  -  lat  (ht  +  Ahf) 


where : 


T 

hi  -  J  C  (T)  dT 
o  "i 


and: 


Cn  *  Zoj  C 
P  1  Pi 


Consideration  will  be  further  restricted  to  the  axis  (x 
and  3/dx  *  0,  reducing  Equations  (1)  -  (5)  to: 


0)  where  u 


|£  +  §£V.0 

at  3y 


(12) 


3pa^  3pva^ 
3t  +  3y 


3ai 

»y 


)  + 


Determination  of  the  chemical  production  term,  w^,  will  employ  the 
method  developed  by  Pratt  This,  in  turn,  uses  the  thermodynamic 
representation  individual  species  developed  by  McBride,  et.al.^  Accordingly, 
a  fifth-order  polynomial  is  used  in  the  present  analysis  for  the  sensible 
enthalpy  and  specific  heat  of  each  species: 


hi 


5 

S  C^  tJ/ 

j-1 


j 


(13) 


Cpi 


5 

i  ct1  tJ_1 

j=i 


(14) 


The  sixth  coefficient  defined  by  McBride,  et.al.^,  is  identified  herein  as  the 
heat  of  formation  at  Ah^.  These  coefficients  are  additionally  used  to  define 
the  equilibrium  constant  and  either  the  forward  or  backward  rate  constant 
(whichever  is  not  otherwise  specified)  from  the  Law  of  Mass  Action. 


Details  of  the  technique  for  determination  of  the  chemical  production  term 
are  provided  in  Reference  14.  For  the  present  discussion  it  may  be  regarded 
as  a  "source"  term  for  Equation  (12)  that  is  a  function  of  the  instantaneous 
pressure,  enthalpy  and  mixture  composition: 

w £  -  f(p,  h,  <Xj)  ;  j  =  1,  NS  and  i  =  1,  NS  (15) 


The  CREK  code  described  in  Reference  14  is  contained  in  the  present 
computermodel  as  an  auxiliary  program  that  integrates  Equation  (15)  by  the 
Newton-Raphson  iteration  technique  described  in  Reference  14: 


^  Pratt,  D.T.,  "Calculation  of  Chemically  Reacting  Flows  with  Complex 
Chemistry"  in  Studies  in  Convection.  B.E.  Launder,  Ed.,  Academic  Press,  NY 
1977. 

^  McBride,  B.J.,  Heims,  S.,  Ehlers,  J.G.  and  Gordon,  S.,  "Thermodynamic 
Properties  to  6000°K  for  210  Substances",  NASA  SP-30001,  1966. 
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t+At 

Aai-reaction  /  —  dt  (16) 

P 

This  integration  is  performed  at  constant  pressure  and  enthalpy,  but  includes 
the  variation  in  temperature,  density  and  composition  due  to  reactions.  (The 
scheme  admits  both  finite-rate  reactions  and  fully-equilibrated  reactions;  only 
the  finite-rate  option  is  of  present  interest.)  Blending  of  the  change  in 
composition  due  to  reactions  with  the  change  due  to  diffusion  will  be  discussed 
later  in  connection  with  the  solution  of  Equation  (12).  However,  it  should  be 
recognized  that  the  changes  in  temperature  and  density  obtained  in  this  step  of 
the  solution  are  incomplete  -  the  complete  changes  are  obtained  from  the 
solution  of  Equation  (10)  and  (11).  The  influence  of  reactions  enters 
implicitly  through  the  use  of  Equations  (6)  and  (7).  If,  for  example.  Equation 
(11)  had  been  written  in  terms  of  temperature  rather  than  enthalpy,  a 
temperature  "source"  term  involving  the  product  w^Ah?  would  appear  and  the 
change  in  temperature  obtained  concurrently  with  the  solution  of  Equation  (16) 
would  be  used  to  represent  time- integration  of  this  term. 


Integral  Form  of  the  Continuity  and  Energy  Equations 


Equations  (10)  and  (11)  can  be  converted  from  a  pair  of  partial 
differential  equations  in  two  independent  variables  (y  and  t)  to  a  pair  of 
ordinary  differential  equations  with  time  (t)  as  the  independent  variable  by 
analytically  integrating  over  y  and  representing  the  integrals  in  terms  of 
shape  parameters  that  are  functions  of  time.  The  limits  of  integration  are  the 
outer  edge  of  the  mixing  layer  on  the  propellant  side,  62,  and  the  outer  edge 
of  the  mixing  layer  on  the  air-side,  6^.  We  choose  to  form  a  pair  of  equations 
by  integrating  first  from  y  =  0  to  y  =  6j  and  again  from  y  =  0  to  y  =  62*  The 
subscripts  0,  1  and  2  will  be  used  hereafter  to  denote  conditions  at  y  =  0,  6^ 
and  62  respectively. 

The  continuity  equation.  Equation  (10),  becomes: 


d 

dt 


o 


pdy  - 


d6i 

dt 


+ 


0  ( i  =  1,2) 


(17) 


where  the  boundary  conditions  v  =  0  at  y  =  0  and  3/3y  =  0  at  y  =  6^  have  been 
imposed  and  L'Hopital's  rule  has  been  observed. 
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The  energy  equation.  Equation  (11),  similarly  becomes: 


phdy  - 


Pivihi 


(f*)  + 

3y  o 


d_ 

dt 


(P«i> 


(18) 


The  velocity  vj  can  be  eliminated  by  multiplying  (17)  by  h*  and  substituting  it 
into  (18),  giving: 


d  .  .  .  d  {i  .  kTo  3h 

Tt{  phdy -hi  jjJ  pdy--^; 


+  5  idt  (i 


1,2) 


(19) 


Evaluation  of  the  integrals  in  Equation  (19)  requires  representation 
of  the  integrands  by  some  appropriate  distribution  functions  satisfying  the 
boundary  conditions  at  y  *  and  $2*  Additional  accuracy  can  be  obtained  by 
also  satisfying  the  differential  equations  ((10  and  (11))  at  y  *  0.  Again, 
continuity  can  be  substituted  into  the  energy  equation  to  eliminate  (3v/3y)Q, 
yielding: 


dho  ,3  ,KT  3hxx  .  d£ 
7Z~  =*  (TT.  (?“  TTJ)  +  77 


°  dt 


1 3y  'cD  3y 


dt 


(20) 


Distribution  Functions 


Let  f(y)  represent  any  of  the  dependent  variables.  The  boundary 
conditions  are: 


f—  *  0  and  f  ■  fd 

3y  1 

at  y  *  5^ 

(21) 

H* 

H 

Hi 

O 

at  y  ■  0 

(22) 

The  following  Fourier  series  is  chosen  to  represent  f(y): 
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f  ■  a  +  b  cos  irq  +  c  cos2Trr) 
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The  primary  role  of  the  distribution  function  is  to  permit 
evaluation  of  the  integral  of  the  function,  viz: 

*1  b(VS2) 

J  fdy  »  a6i - - -  sin  irn0 - ^ — L  sin  2irn0 


and  the  derivative  of  the  integral,  viz: 


'll 

J  fdy  m  a^  +  Sja  + 


Evaluation  of  (37)  is  aided  by  the  following  definitions: 
Yi  (a)  -  a 


Y2(a,  b,  c.  Si)  *  a6i  - 


in  irn0  "  — 27 


Y3(a2,  a3 ,  6)  *  -(«!  -  52)(a2cosirn0  +  a3cos2rn0)5 
I  (a,  b,  c)  ■  -(bsinvnQ  ■  i  c  sin  2vn0)/Tr 


where 


(«!  -  «2)2 


i  -  1 


+  - - -  i  -  2 

(«i  -  *2)2 

Using  these  functions.  Equation  (37)  becomes: 


J  fdy  ■  Yi(a)6i  +  y 2(a,  b,  c,  6^  +  I  (a,  b,  c)(6x  -  62) 


+  Yj(b»  c,- 
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where  a,  b,  c  are  given  by  Equations  (33),  (34)  and  (35),  and 
Y2(a,  b,  c,  62)  *  Y2(ai*  bi*  ci»  ^  • 

+  Y2(a5»  b5,  c5,  5i)6J 
Accordingly,  (43)  can  be  written  as: 

6i 

^  j  fdy  *  [Yi(a)]6i  +  [Y2(ai.  bi»  ci*  5i))f0 
o 

^  [y2(a2®  P2*  c2,  5 3 f i 

[Y2(a3.  bj,  C3,  6 ^ ) ) f 2 

+  [I(a,  b,  c)  +  Y2(a4»  b4,  c4,  6^  +  X3(b,  c, 


(44) 


(45) 


(«!  -  62)2 


2)16X 


+  [-I(a,  b,  c)  +  Y2(a5»  bg,  C5,  6^)  +  X3(b,  c. 


-61 


(61  -  62): 


-2)]65 


Finally,  the  distribution  function  is  also  used  to  evaluate  the  first  and 
second  derivatives  with  respect  to  y  at  y  =  0: 


(f?l*  *7^,  <‘ !l"  •%  * 2 ' si"  ”„) 

-IT2  /,  .  „  \ 

-  »  -  (b  cos  nn0  +  2  c  cos  nn0) 

Uy2Jo  _  6z>2 


(46) 

(47) 


Species  Equations 

The  sane  integration  steps  as  performed  above  could  be  applied  to  the 
species  equations.  Equation  (12), which  would  yield  a  system  of  equations 
similar  to  Equation  (19)  but  involving  a  separate  limit  of  integration  for 
the  two  layers  (i  ■  1,2)  and  each  species  (j  *  1,  NS).—  However,  consistent 
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with  the  assumption  of  unit  Prandtl  and  Lewis  numbers,  it  may  be  assumed  that 
each  species  will  diffuse  to  the  "edge"  of  the  mixing  layer  as  defined  by  the 
conduction  of  heat  across  the  layer,  i.e.,  by  the  thermal  thickness  6^. 
Equating  these  thicknesses  makes  the  integral  form  of  the  species  equations 
redundant,  and  it  becomes  necessary  to  solve  the  species  equation  only  at  the 
interface  location  y  =  0.  Each  of  the  species  mass  fractions  will  be  described 
by  a  distribution  function  f^  of  the  form  given  by  Equation  (23),  viz: 


^i  =  ai  +  ^i  cos  UT1  +  ci  cos  2irn  (i  =  1*  NS) 


where: 


ai  “  ai  K0’ail’ai2*5i*  52>  .  (49> 

bt  =  bt  (ailfoi2)  (50) 

ci  =  ct  (aio,ail,ai2,61,62)  (51) 

and  the  coefficients  a^,  b^,  Cj  have  the  same  properties  as  displayed  in 

Equations  (30)  -  (35)  and  (46),  (47). 

The  species  mass  fractions  at  the  interface  are  determined  from  Equation 
(12)  evaluated  at  y  =  0: 


dai0  r-3  f  _  3 ai>,  .  . 

Po  —  "Li*  [  pDT3TJJ0+ 


Final  System  of  Enerev  and  Soecies  Eauation 


Equations  (19),  (20)  and  (52),  aided  by  (45),  (46)  and  (47),  form  a 
coupled  linear  system  of  first-order  ordinary  differential  equations.  The 
coupling  occurs  through  the  equations  of  state,  Equation  (6)  and  (7).  In 
particular,  the  enthalpy  at  the  interface,  hQ,  yields  the  temperature  at  the 
interface  by  iterative  solution  of  (7).  The  coupling  can  be  simplified,  and  the 
iteration  eliminated  by  differentiating  the  equations  of  state  and  using  the 
density,  rather  than  enthalpy,  as  the  dependent  variable  for  the  energy 
equation: 


GASL  TR-276 


P/p  -  p/p  +  T/T  +  MW  (Z  (d^MWi) 
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h  -  £  a^hi  +  Ah?)  +^ai(E  C^TJ*1)! 
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where 


MW  =  [Zc^/MWi)]' 
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Combining  (20),  (53)  and  (54)  yields: 
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Equation  (56)  now  replaces  (20),  and  the  system  (19),  (52)  and  (56) 
can  be  expressed  in  matrix  form  as: 
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where: 

i  =  1,  NS  +  3  (number  of  equations) 
j  «  1,  NS  (number  of  species) 

Ax,i  ■  0 
A1>2  =  0 
Ai t 3  *  C 

A,.,  +  j  -  P»[  l  Ah?  *\  (l  - 

"  fjr  — ^  +  (1  -  Po*.,.)  P  -  Si4  5 
^  p  3y2Jo  i 

^2  ^  i  =  ^3)  A2(a^4, 

+  A 3 ( a2 >  3 

A2  >  2  *  ’^(^1  >  ^2>  ^3)  ^  ^2^  »  £3^ 

+  A3(a2 »  ^3 1  ~6i/( 

A  2  >  3  =  ^2^»  ^  3  *  »  ^  1  )  ( A 1  >  3  ) 


A2,3+j  *  (£d  ‘  Siai)  CijTJ)A2(a1l,  a2* 

-  h1A2(a1,  b1,  c1,  5j) 
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B2  *  ^2^al^»  a2^’  a31*  ^1^  (^Al,3  +  ^1^  P  (71) 

'(S  3y  >0  -Plhlx2(al2*  a22*  a32»  61>-*- 
-P2h2X2(a13,  a23,  a33,  Si) 

A3>1  =  I  (alt  a2,  a3)  +  X^a^,  a24»  a34,  62)  +  ...  (72) 

A2  ^  3  *  X^(sj)  ~  X  (ax ,  a2 1  a3)  "f*  X2(a2  .  a2  9  a3  >  62)  •••  (73) 

A3,3  =  X2(a11,  a2 1 t  aj1,  S2)  (hQ  +  A1>3)  +  (74) 


A3,3+j  *  (Jd  '  Cijlj)X2(a11,  a^,  a31,  S2 ) 

-  h2X2(a1,  b1,  c1,  62) 

B3  *  X2(a11,  a^,  a31,  62)  (^A1>3  +  62)  p 

~(^  |y  )Q  -(Ph)lX2(ax2,  a22,  a32,  «2)... 
-(ph)2X2(a13,  a23,  a33,  62).  .  . 

and: 

A3+j,l  m  0 
a3+J,2  *  0 
a3+J,3  *  0 

A3+j , , 3+m  *  6jm  (Kronecker  delta) 
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(75) 


(76) 


(77) 

(78) 

(79) 

(80) 
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where  j  =  1,  NS  and  m  =  1,NS. 

In  view  of  the  fact  that  A3+j  ,  3^  =  6jm  (i.e.,  0  for  j  *  m  and  1  for  j  = 
m)  it  is  evident  that  the  species  equations  are  uncoupled  and  can  be  integrated 
independently  of  the  remainder  of  the  system  (over  a  single  time  step): 


a,  =  B 


i  "  fl3+i 


In  this  case  the  B^'s  can  be  redefined  as: 


B'i  =  Bi  -  I  ^i,3+j  aj  i  1>2,3 

Furthermore,  since  Aj_j  =  A32  =  Os 
Po  =  B'l/Al,3 

The  B*2  and  B*  3  can  be  redefined  again  as: 

b"2  =  b'2  -  a2>3p0 


B  3  =  B  3  -  A3>3p0 
and  the  remaining  2  by  2  system  is: 


a2,1  51  +  a2,262  =  b"2 


a3, 1  *1  +  a3,2°2  "  B  3 
which  is  easily  solved: 


•  B"2A3  2  -  B"3a2  2 

1  a2, 1a3,2  '  a3 , 1A2, 2 


m 

m 

m, 

V, 
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5 


2 


b<I3A2.1  -  b<2A3 , 1 
a2, 1A3,2  '  a3 , 1a2 , 2 


(90) 


Thus  the  final  system  of  differential  equations  consists  of  (81),  (83), 
(87)  and  (88).  These  are  integrated  over  a  time  step  by  a  standard  Runge-Kutta 
procedure. 


Turbulence  Modeling 

Three  options  are  provided  in  the  computer  code  for  representation  of  the 
thermal  conductivity  k^  and  species  diffusivity  Dip.  It  should  be  noted 
that  kT  only  appears  in  ratio  to  Cp  and  D-p  only  appears  in  a  product  with  p. 
Hence,  for  unit  Prandtl  and  Lewis  numbers: 

^/Cp  *  pD-p  *  p  (91) 

The  first  option  is  to  use  the  laminar  viscosity  provided  by 
Sutherland's  law: 


U  =  Sx  W2/  (T  +  Sj)  (92) 

where  Sp  and  S2  are  constants  in  the  appropriate  set  of  units  (e.g.,  Sp  *  2.27 
x  10”®  lb  sec/ft^/°R2  and  S2  *  198. 6°R). 

The  second  option  is  to  use  a  modified  form  of  eddy  viscosity;  the 
following  is  proposed  here  (based  on  analogy  with  accepted  kinematic  models): 


P  ■  S3  W  hp'  6p  (pp  -  pg)  +  $2  (p2  "  P0»  (93) 

where  S3  is  an  empirical  constant.  (We  have  used  S3  *  .0143  lb  sec/ ft^  in  the 
code . ) 

The  third  option  is  to  use  a  model  of  the  two-equation  turbulent 
thermal  energy  equation.  The  model  is  derived  by  analogy  with  the  widely-used 
two-equation  model  of  the  turbulent  kinetic  energy  equation.  The  two  variables 

in  the  present  model  are  the  mean  square  temperature  fluctuation  82  and  a 
thermal  dissipation  coefficient,  <p.  The  thermal  conductivity  is  related  to 
these  variables  by: 
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j 

&T  "  copCp02/$ 


(94) 


These  two  variables  are  governed  by  differential  equations  developed  from  the 
turbulent  thermal  energy  equation  by  analogy  with  the  model  equations  for  the 
turbulent  kinetic  energy  equation: 


D62  kT  ,3Tx2 
p  Dt  “  Cp  (3y' 


3_ 

3y 


30£v 

^Cp  3y  ; 


+  p+ 


0 


Dt 


_  3T 

lS  7  * 


c2 


3_  ,^T  3£ 
3y  *Cp  3y 


) 


(95) 


(96) 


We  have  assumed  that  both  6^  and  cp  can  be  described  by  the  same 
distribution  function  used  for  the  other  variables.  Equation  (23), which, 
however,  implies  that  their  ratio  is  independent  of  y,  and  hence  k<p/Cp  and  pD-p 
are  likewise  independent  of  y.  It  is  therefore  only  necessary  to  evaluate 
Equation  (93)  and  (94)  at  the  interface,  y  »  0,  to  obtain  the  differential 
equations  for  6-QJ  and  $q: 


(97) 


(98) 


In  this  case.  Equations  (95)  and  (96)  are  added  to  the  system  of 

equations  to  be  integrated,  and  0Q2  and  $q  are  added  to  the  list  of  dependent 
variables. 

Chemical  Kinetics  Modeling 


The  final  form  of  the  species  equation, Equation  (82),  can  be  written 
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a  =  a  +  a 

i  i  i 

diffusion  reaction  (99) 


where  a.  and  a.  are  identifiable  with  the  two  terms  in 

diffusion  .  ^-reaction 

B3+i»  Equation  (81).  Equation  (97)  can  be  integrated  in  a  variety  of  ways, 
offering  various  advantages  in  speed,  simplicity,  stability,  etc.  We  have 
found  the  most  stable  procedure  to  be  to  evaluate  dj  first,  by  calling 

reaction 

the  CREK  code  to  integrate  the  chemical  production  rate  over  the  time  step  At. 
Equation  (82)  is  then  integrated  by  the  same  Runge-Kutta  scheme  as  the 
continuity  and  energy  equations,  yielding  the  combined  effects  of  diffusion  and 
reaction.  The  permissible  step  size  At  is  governed  by  an  error  bound  applied 
to  all  equations  integrated  by  the  Runge-Kutta  procedure.  However,  when  the 
rate  of  chemical  reaction  becomes  very  fast,  the  integration  procedure  in  CREK 
can  be  stable  but  yield  changes  in  composition  that  are  too  large  for  the 
Runge-Kutta  procedure  to  integrate  accurately.  Therefore,  no  more  than  a  180°R 
change  in  temperature  due  to  reaction  is  permitted  in  any  time  step.  If  this 
limit  is  exceeded  the  time  step  is  halved  and  the  evaluation  of  reaction  ’•s 
repeated  before  proceeding  to  the  Runge-Kutta  integration. 

The  CREK  code  is  very  versatile,  as  well  as  accurate  and  stable.  The 
chemical  constituents  of  the  mixture  are  specified  by  their  chemical  names,  the 
polynomial  coefficients  defining  their  enthalpy,  Cjj,  and  the  atomic  weights  of 
the  elements  making  up  the  constituents.  The  kinetic  rate  mechanism  is  defined 
symbolically  by  the  names  of  the  reactants  and  products  in  each  reaction  and 
either  the  forward  or  backward  rate  constant  in  the  form 


R 


kt  *  10AiTBiexp  (-Ci/T)  (100) 

All  this  data  is  specified  in  an  input  file  and  is  therefore  easily  modified 
(without  any  recoding). 
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4.0 


DESCRIPTION  OF  THE  COMPUTER  CODE 


Overview 


The  preceding  analysis,  while  complicated  in  derivation,  results  in  a 
fairly  compact  and  efficient  computer  code.  The  Fortran  code  UMLAC  (Unsteady 
Mixing  Layer  and  Chemistry)  executes  on  an  HP  1000  minicomputer  using  the  EMA 
(Extended  Memory  Area)  capability  of  this  machine.  Although  the  calls  to  the 
CREK  program  are  inefficient  (the  main  program  transfers  data  to  CREK,  then 
suspends  itself  while  CREK  executes,  (including  a  complete  initialization  of 
the  thermodynamic  and  kinetic  rate  data  for  each  call,)  CREK  passes  its  results 
back  to  the  main  program,  which  then  resumes  execution,)  the  running  times  are 
not  unreasonable  (i.e.,  10's  of  minutes,  typically.)  A  modest  amount  of 
recoding  to  include  the  CREK  code  as  a  true  set  of  sub-routines  and  eliminate 
the  reinitialization  would  substantially  accelerate  the  execution  speed. 

Input  data  is  contained  on  three  formatted  ASCII  files;  one  is  read  by  the 
main  program  on  logical  unit  5,  the  second  is  read  by  the  main  program  on 
logical  unit  9,  and  the  third  is  read  by  both  the  main  program  and  by  CREK  on 
logical  unit  15.  The  first  file  contains  the  initial  conditions  and 
integration  parameters  for  the  case  at  hand.  The  second  contains  the  boundary 
condition  data,  i.e.,  pressure,  inviscid  flow  densities  on  the  air-side  and 
propellant-side  of  the  interface,  and  species  mass  fractions  on  both  sides,  all 
as  functions  of  time.  The  third  contains  the  thermodynamic  and  kinetic  rate 
data;  the  thermodynamic  data  is  read  by  both  UMLAC  and  CREK. 

In  view  of  the  importance  of  the  inviscid  entropy  layer  on  the  propellant- 
side  of  the  interface  pointed  out  in  the  earlier  discussion,  the  program 
permits  up  to  11  streamlines  in  the  propellant-side  flowfield  to  be  specified 
as  part  of  the  boundary  condition  data.  The  local  "edge"  properties  are 
obtained  by  linear  interpolation  in  the  y  direction  using  the  value  of  62  at 
the  current  time  step,  as  indicated  in  Figure  (2).  Since  the  task  of  preparing 
this  massive  amount  of  boundary  data  is  quite  formidable,  a  pre-processor  code 
DCREK  has  been  written  to  provide  an  automated  link  between  the  inviscid  muzzle 
blast  flowfield  code  DAWNA*®  and  the  present  mixing  layer  code  UMLAC.  DCREK 
reads  the  pressure  history  along  the  interface  from  DAWNA  and  the  densities 
along  each  of  the  selected  streamlines.  (The  density  is  isentropically 
adjusted  to  the  interface  pressure  should  a  pressure  difference  exist  between 
the  streamline  and  the  interface).  It  also  reads  the  initial  chemical 
composition  for  the  streamlines,  which  can  be  equated  to  the  equilibrium 
composition  of  the  propellant  gases  at  the  muzzle  exit,  since  that  composition 
remains  effectively  frozen  through  the  subsequent  supersonic  expansion  and 


*®  Ranlet,  J.  and  Erdos,  J.  "Description  of  Fortran  Program  DAWNA  for 
Analysis  of  Muzzle  Blast  Flowfield"  BRLCR-302,  March  1975. 
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recompression  across  the  Mach  disc.  DCREK  converts  the  pressure  and  density  to 
an  enthalpy  history  for  each  streamline  using  the  initial  composition,  and  then 
performs  a  finite-rate  calculation  along  the  streamline  (preserving  the 
pressure-enthalpy  history).  DCREK  stores  the  pressure,  composition  and  (post¬ 
reaction)  density  history  for  each  streamline  on  an  output  file  that  becomes 
the  input  data  file  for  UMLAC.  The  same  calculation  is  also  performed  on  the 
air-side,  but  only  for  the  interface  streamline.  DCREK  automatically 
initializes  the  air-side  composition  to  a  mixture  of  76.55%  N2  and  23.45%  O2 
(by  mass). 

Since  the  same  thermodynamic  and  kinetic-rate  data  file  is  used  by  CREK  in 
both  the  DCREK  and  UMLAC  codes,  it  should  include  the  appropriate  air 
constituents  and  air  reaction  mechanisms  as  well  as  the  propellant  gas 
constituents  and  mechanisms.  For  typical  muzzle  blast  conditions  the  air 
constituents  are  N2,  N,  O2,  0  and  NO,  and  the  reactions  among  them  are  oxygen 
dissociation  and  the  Zeldovich  mechanism  for  nitric  oxide  formation: 


02  +  M 

t 

0  +  0  +  M 

(101) 

n2  +  0 

t 

NO  +  N 

(102) 

02  +  N 

t 

NO  +  0 

(103) 

Input  data 


The  initial  conditions  and  option  control  data  for  UMLAC  are  given  in 
4  +  NS  number  of  lines,  where  NS  is  the  number  of  species.  The  format  and 
variables  are  shown  in  Table  1,  and  the  variables  are  defined  as  follows: 

Line  1:  Cl  Initial  value  of  time  step.  At 

ER  Max  integration  error  in  any  variable  per  time  step 

FCT  Initial  value  ratio  6^/(6^  +  62) 

Line  2:  TINT  Initial  value  of  time 

DELINT  Initial  value  of  total  thickness,  +  62 

TFINAL  Final  time  (at  which  run  is  terminated) 

DINF  Reference  length  (e.g.,  gun  barrel  diameter) 

THETIN  Initial  value  of  0^  (if  used) 

DHINT  Initial  value  of  9  (if  used) 

Initial  value  of  time  step-counter  (e.g.,  0) 


Line  3: 


KINT 
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KFINAL 

LINT 

ITURB 


Maximum  value  of  time  step  counter 
Print  interval  counter 

Turbulence  model  flag:  -1  for  the  two-equation  TTE 
model,  0  for  the  eddy  viscosity  (conductivity)  model, 
+1  for  laminar  viscosity  (Sutherland's  law) 


ICHEM 


IWRITE 


Chemistry  flag:  0  for  frozen  flow,  1  for  finite-rate 

kinetics,  2  for  fully  equilibrated  chemistry. 

Debug  flag:  0  for  normal  printout,  1  for 

comprehensive  data  dumps 


Line  4:  IALP 

IELE 


Number  of  species,  NS 
Number  of  elements 


Lines  4  +  NS: 

MW(I) 

RNAME(I) 

DUM 


Molecular  weight  of  i^  species 
Chemical  name  of  ifc®  species 
Dummy  variable;  not  used 


In  addition,  certain  semi -permanent  data  is  stored  in  DATA  Statements 
contained  in  a  BLOCK  DATA  subroutine.  The  variables  and  their  currently 
assigned  values  are: 


CO,  Cl,  C2,  C3 


0.09,  1.43,  0.77,  1.92 


2117.  psf 


RHOINF 


0.002376  slugs/ft3 


GAMINF 


CPINF 


6006.  ft2/sec2/°R 


The  nondimensional  constants  are  the  values  of  Cq,  Cj,  C2  and  C3  appearing 
in  the  two-equation  TTE  turbulence  model.  The  values  of  pressure, density, 
ratio  of  specific  heats  and  specific  heat  at  constant  pressure  correspond  to 
standard  atmospheric  conditions.  Certain  variables  transferred  into  UMLAC  from 
DCREK  may  be  non-dimens ionalized  with  respect  to  these  parameters,  as  well  as 
the  reference  length  DINF  appearing  in  Line  2  of  the  input  data  file.  The 
values  of  these  parameters  in  entirely  arbitrary,  and  they  may  be  all  defined 
as  1.0,  for  example,  or  may  be  defined  in  S.I  units,  etc.  The  only  important 
consideration  relative  to  their  definition  is  that  they  are  assigned  consistent 
values  in  DCREK  and  UMLAC.  (In  this  respect,  values  of  1.0  may  be  the  best 
choice. ) 


r  laVl 


The  boundary  conditions  are  given  on  a  second  input  file,  as  described  in 
Table  2.  In  general,  this  file  will  be  automatically  generated  by 
Program  DCREK  and  require  no  alterations  by  the  user.  However,  it  can  be 
manually  prepared,  if  desired.  The  variables  are  defined  as  follows: 


Line  1:  I ALP  Number  of  Species,  NS 

NSL  Number  of  streamlines  on  side  #2  (max  of  11) 

ISCT(l)  Number  of  dummy  entries  in  the  data  for  the  i^ 

streamline  (corresponds  to  the  time  it  takes 
before  the  ifc^  streamline  crosses  the  Mach  disc). 

Line  2:  TFIT  Time  (at  which  the  following  values  occur) 

P  Pressure 

RH02(1)  Density  on  side  ill  of  the  interface 

RH01  Density  on  side  II 1  of  the  interface 

RH02(I)  Density  on  ifc^  streamline  on  side  111  (I“2,  NSL) 

Line  3:  ALPHA2(1,I)  Mass  fractions  of  species  #1  through  NS  on  side 

111  of  the  interface. 

Line  4:  ALPHAl(I)  Mass  fractions  on  side  #1  of  the  interface. 

Line  5:  ALPHA2(J,I)  Mass  fractions  on  streamline  on  side  II 2. 

***  The  order  of  the  species  must  match  that  used  in  lines  5  through  4  + 
NS  on  the  input  data  file  on  logical  unit  5  (Table  1)  and  that  used 
in  the  THERMO  section  of  the  CREK  data  file  (Table  3).  *** 

Line  6:  ZSL(J)  Axial  distance  from  interface  to  the 

streamline. 


The  thermodynamic  and  reaction  kinetic  data  is  on  a  third  input  file, 
which  is  described  in  Table  3.  This  file  is  read  by  the  Program  CREK  which  is 
called  by  UMLAC,  and  the  identical  file  should  be  used  by  DCREK  to  generate  the 
boundary  condition  data.  The  variables  in  this  file  are  defined  as  follows: 


Line 

1: 

ELEMENTS 

Key  word  identifying  the  type 
follows 

of 

data 

which 

Line 

2: 

NAME  (I) 
ATWT(I) 

VAL  (I) 

Chemical  name  of  the  it*'  element 
Atomic  weight  of  the  it^'  element 
Valence  of  the  i^  element 

Line 

3: 

THERMO 

Key  word  identifying  the  type 
follows  — 

of 

data 

which 
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Line  A: 


Line  5: 


NAME  (I) 
T1(I) 

T2(I) 


Chemical  name  of  i^1  species 

Upper  limit  of  temperature  range  for  first  set  of 
polynomial  coefficients 

Upper  limit  of  temperature  range  for  second  set 
of  polynomial  coefficients 


C(I,J)  7  polynomial  coefficients,  as  defined  by 

Reference  (10),  for  the  first  temperature  range, 
followed  by  7  more  for  the  second  temperature 
range. 


Line  6:  MECHANISM  Key  work  identifying  type  of  data  which  follows 


Line  7s 


D,E,F,G,H,I  Chemical  names  of  reactants  and  products.  M 
denotes  a  third-body  which  can  be  any  species. 


A,B,C 


Terms  in  the  rate  constant  given  by  k“10^T®  exp 


(-C/T) 

K  K  can  be  either  a  blank  or  one  of  the  following  A 

character  expressions:  REVE  denotes  that  reverse 
rate  data  has  been  specified,  CGS  denotes  that 
the  data  is  not  in  S.I.  units,  COMM  denotes  a 
comment  line  (not  interpreted  as  data).  Any 
other  expression  will  be  ignored. 

Sample  data  sets  for  the  input  data,  boundary  condition  data,  and  thermo 
and  kinetic  data  files  are  presented  in  Figures  (3),  (A)  and  (5). 


Output  Data 


The  output  from  UMLAC  is  divided  into  three  categories.  The  first  is  a 
repetition  of  input  data  (in  card  image  format)  to  assist  in  checking  for  input 
errors.  The  record  is  a  set  of  statements  of  calculated  constants  and  initial 
conditions.  The  third,  which  comprises  most  of  the  output,  is  a  statement  of 
the  solution  and  boundary  conditions.  A  sample  output  is  provided  in  Figure  6. 

The  first  block  of  data  is  the  repetition  of  the  data  read  from  LU  'It  5, 
and  matches  the  description  given  in  Table  1.  The  initial  part  of  the  output 
data  shown  in  Figure  6  matches  the  input  data  shown  in  Figure  3. 

Following  this  is  a  short  summary  of  calculated  scale  factors  used  to 
nondimensionalize  certain  variables  in  the  program,  and  a  message  indicating 
that  the  data  file  on  LU  II  9  (Figure  A)  has  been  successfully  read.  Thereafter 
is  a  synopsis  of  the  streamline  data  showing  the  initial  time  and  postion  at 
which  each  streamline  begins.  The  warning  statement  at  the  end  of  this  set  of 
data  cautions  the  user  that  if  the  program  continues  beyond  0.02898  seconds  it 
will  use  the  boundary  data  at  that  time  to  the  final  time,  0.10  seconds  in  this 
example.  — 
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The  set  of  calculated  initial  conditions  for  the  species  are  given  next, 
in  the  same  format  as  used  previously  for  the  species  input  data  on  LU  if  5.  In 
this  instance  the  first  two  columns  of  data  now  contain  the  mass  fractions  of 
the  species  on  the  air  side  of  the  interface  (i.e.  at  y  *  5^)  and  on  the 
propellant  side  (i.e.  at  y  »  6t).  The  latter  values  have  been  interpolated 

from  the  streamline  data  at  the  initial  time. 

The  solution  data  begins  with  a  set  of  labels  identifying  the  current 
value  of  the  time  set  counter,  K.  The  following  glossary  identifies  the  output 
variables  in  the  order  they  appear  : 

K  Time  step  counter 

TIME  Current  time  (seconds) 

DELTA  1  Thickness  of  air-side  layer;  6^  (ft) 

DELTA  2  Thickness  of  propellant-side  layer;  62  (ft) 

THETA  **  2  Mean-squared  temperature  fluctuation;  0^ 

PHI  Thermal  energy  dissipation  parameter;  $ 

RHO  Density  at  y  *  0;po  (non  dimensional) 

KTO  Thermal  conductivity  at  y  •*  0 

P  Pressure  (non  dimensional) 


TEMPI 


TEMP2 


YZERO 


ALPHA 


Density  at  y  =  6^; 

Density  at  y  =  62;  P2 
Temperature  at  y  ■  6j,  Tj  (°R) 

Temperature  at  y  =  62*  ^2  (°R) 

Location  of  62  in  streamline  coordinate  system 
Temperature  at  y  *  0,  TQ  (°R) 

Enthalpy  at  y  *  0,  hQ  (ft2/sec2) 

Array  of  species  mass  fractions  at  y  *  0; 

(same  order  as  given  by  the  user  in  the  input  data 
files. ) 

dp/dt 
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DRH1DT 

(dph/dt)  at  y  = 

DRH2DT 

(dph/dt)  at  y  *  $2 

ALPHA 1 

Mass  fractions  at  y  *  6^; 

DALP1DT 

(da^/dt)  at  y  a  ij 

ALPHA2 

Mass  fractions  at  y  *  62;  ai2 

DALP2DT 

(da^/ dt)  at  y  *  62 

Following  each  solution  statement  is  a  set  of  3  integer  numbers  that 
relate  to  the  integration  algorithm.  The  normal  values  are  4,  1  and  1.  Also 
shown  here  are  the  current  values  of  the  (nondimensional)  time  step  and  the 
values  that  will  be  attempted  on  the  next  step. 

The  solution  data  will  be  repeated  until  either  K  >  KFINAL  or  TIME  > 
TFINAL . 

The  presence  of  muzzle  flash  is  characterized  by  the  occurrence  of  a  rapid 
increase  in  the  interface  temperature  (TEMP),  which  otherwise  drops  off  as  the 
pressure  (P)  decays  and  tends  to  remain  bracketed  by  the  bounding  temperatures 
(TEMPI  and  TEMP2).  The  species  composition  at  the  interface  (ALPHA)  provides 
additional  diagnostic  information  regarding  the  chemical  reactions  leading  to 
or  suppressing  flash.  In  general,  the  user  is  advised  to  examine  the  time 
variations  of  pressure  (P)  and  temperature  (TEMP)  to  determine  if  flash  occurs 
in  a  particular  case,  and  then  examine  the  time  histories  of  the  species 
(ALPHA),  the  thicknesses  of  the  two  layers  (DELTA1  and  DELTA2),  and  the 
bounding  temperatures  (TEMPI  and  TEMP2)  for  diagnostic  information. 
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5.0  SUMMARY 

A  model  of  one-dimensional,  unsteady  turbulent  mixing  and  chemical 
reaction  along  the  interface  between  the  propellant  gas  expelled  from  a  gun 
muzzle  and  the  blast-heated  air  has  been  developed.  The  model  incorporates  all 
the  essential  aspects  of  the  unsteady  flowfield  processes,  but  the 
computational  burden  is  economized  considerably  by  employing  a  two-layer 
integral-method  for  solving  the  governing  gas -dynamic  equations.  Derivation  of 
the  integral -method  of  analysis  and  of  the  resulting  system  of  coupled,  first 
order  ordinary  differential  equations  has  been  presented  in  detail.  Assumption 
of  unit  Lewis  and  Prandtl  numbers  and  of  a  single  definition  of  the  "edge"  of 
the  mixing  layer  for  enthalpy  and  for  species  diffusion  nearly  uncouples  the 
energy  equation  from  the  species  conservation  equation,  and  requires  that  the 
species  equation  be  solved  only  at  the  interface  (at  no  small  savings  in 
computer  time). 

The  general  requirements  for  operating  the  computer  code  and  the  specific 
input  and  output  data  formats  and  a  glossary  of  Fortran  variables  have  been 
presented. 
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TABLE  1 


DATA  FORMAT  FOR  UMLAC  INPUT  FILE 
(Read  on  Logical  Unit  5) 


4  +  IALP 


Format 

Variables 

3E  10. 

Cl,  ER,  FCT 

6E10. 

TINT,  DELINT 

615 

KINT,  KFINAL 

215 

IALP,  IELE 

3E10..A4, 

DUM,  DUM,  MW 

5X.A4 

II 

II  II  1 

II 

II  II  1 

TABLE  2 


DATA  FORMAT  FOR  UMLAC  BOUNDARY  CONDITION  FILE 
(Read  on  Logical  Unit  9) 


Format 


Variables 


1  1215  IALP,  NSL,  (ISCT(I) ,  1=2,  NSL) 

***  1  +  (IALP-l)/6  lines  are  skipped  *** 

2  6E12  TFIT,  P,  RH02(l),  RH01,  (RH02(I),  1=2,  NSL) 

3  6E12  (ALPHA2( 1,1),  I  *  IALP) 

4  6E12  ALPHAl(I),  1*1,  IALP 

5  6E12  (ALPHA2(J,I),  1=1, IALP) 

***  Line  5  is  repeated  for  J*2,  NSL  *** 

6  6E12  (ZSL(I),  1=1, NSL) 

***  Lines  2  through  6  are  repeated  up  to  91  times.  The  READ  loop  will 
terminate  upon  encountering  an  end-of-file  mark.  *** 
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TABLE  3 


DATA  FORMAT  FOR  CREK  INPUT  FILE 
(Read  on  logical  unit  15) 

Line  Format  Variables 

1  12A4  ELEMENTS,  a  key  name 

2  A2,  7X,  2F10  NAME  (I),  ATWT(I),  VAL(I) 

***  Line  2  is  repeated  for  each  element  in  the  chemical  system.  A  blank 

line  signals  the  end  of  the  elemental  data  *** 

3  12A4  THERMO,  a  key  name 

4  A12,  27X,  2F10  NAME  (I),  T1(I),  T2(I) 

5  5E15  (Cx  (I,J),  J-1,5) 

6  5E15  (Cx  (I,J),  J-6,7)  (C2(I,J),  J=l,3) 

7  5E15  (C2  (I.J),  J=4 , 7 ) 

***  Lines  4-7  are  repeated  for  each  species  in  the  chemical  system.  A 
blank  line  signals  the  end  of  the  chemical  data.  Note  that  the  number,  names 
and  order  of  the  species  must  agree  with  the  data  specified  in  other  input 
files.  *** 

8  12A4  MECHANISM,  a  key  name 

9  12A4,  3F8.3,  2A4  D,E,F,G,H,I,A,B,C,K 

***  Line  9  is  repeated  for  each  reaction.  It  can  be  used  to  specify  either 
the  forward  or  backward  rate  of  a  reaction,  in  which  case  the  opposite  rate  is 
calculated  from  the  equilibrium  constant,  or  it  can  be  used  to  specify  both 
forward  and  backward  rates  (2  lines  per  reaction).  A  blank  line  signals  the 
end  of  the  kinetics  data.  *** 
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l.E-5  1.0E-3  0.5 

.00052  .0200  0.1  .5085  2.704E+01  1.020E+08 


0  9000  01  -0 

1  0 

11  5 
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2345E-00  . 1355E-08 
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39.1000 
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56.0100 

KOH 

FIGURE  3.  TYPICAL  INPUT  DATA  FILE  FOR  UMLAC  (LOGICAL  UNIT  5) 
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ELEMENTS 

C  12.01115  4.000000 

H  1.007970  1.000000 

0  15.999400-2.000000 

N  14.006700  3.0 

K  39.100000  1.0 

THERMO 

CO  J  9/65C  1.0  1.00  0.00  O.G  300.000  5000.000 

0.29840689E  01  0.14891387E-02-0.57899678E-06  0. 10364576E-09-0.69353499E-14 
-0. 14245227E  05  0.63479147E  01  0.37100916E  01-0. 16190964E-02  0.36923584E-05 
-0.20319673E-08  0.23953344E-12-0. 14356309E  05  0.295S5340E  01 
C02  J  9/65C  1.0  2.00  0.00  O.G  300.000  5000.000 

0.44608040E  01  0.30981717E-02-0. 12392566E-05  0.22741323E-09-0. 15525948E-13 
-0.48961438E  05-0.98635978E  00  0.24007788E  01  0.87350905E-02-0.66070861E-05 
0.20021860E-08  0.63274039E-15-0.48377520E  05  0.96951447E  01 
H  J  9/65H  1.00  0.00  0.00  O.G  300.000  5000.000 

0.25000000E  01  0.0  0.0  0.0  0.0 

0.25471625E  05-0.46011758E  00  0.25000000E  01  0.0  0.0 

0.0  0.0  0.2547162SE  05-0.46011758E  00 

H2  J  3/61H  2.0  0.0  0.0  O.G  300.000  5000.000 

0.31001883E  01  0.5U19458E-03  0.52644204E-07-0.34909964E-10  0.36945341E-14 
-0.87738013E  03-0. 19629412E  01  0.30574446E  01  0.26765198E-02-0.58099149E-05 
0.55210343E-08-0. 181 22726E- 11-0. 98890430E  03-0.22997046E  01 
H20  J  3/61H  2.0  1.00  0.00  O.G  300.000  5000.000 

0.27167616E  01  0.29451370E-02-0.80224368E-06  0. 10226681E-09-0.48472104E-14 
-0.29905820E  05  0.66305666E  01  0.40701275E  01-0.11084499E-02  0.41521180E-05 
-0. 29637 404E-08  0.80702101E-12-0.30279719E  05-0.32270038E  00 
N2  J  9/65N  2.0  0.0  0.0  O.G  300.000  5000.000 

0.28963194E  01  0.15154863E-02-0.57235275E-06  0.99807385E-10-0.65223536E-14 
-0.90586182E  03  0.61615143E  01  0.36748257E  01-0. 12081496E-02  0.23240100E-05 
-0 . 63217520E-09-0 . 22577253E- 12-0 .1061 1587E  04  0.23580418E  01 
0  J  6/620  1.00  0.00  0.00  O.G  300.000  5000.000 

0.25420580E  01-0.27550603E-04-0.31028029E-08  0.45510670E-11-0.43680494E-15 
0.29230801E  05  0.49203072E  01  0.29464283E  01-0.16381664E-02  0.24210303E-05 
-0. 16028432E-08  0.38906964E-12  0.29147641E  05  0.29639931E  01 
OH  J  3/660  l.H  1.00  0.00  O.G  300.000  5000.000 

0.29106417E  01  0.95931627E-03-0. 19441700E-06  0.13756646E-10  0.14224542E-15 
0.393538UE  04  0.54423428E  01  0.38375931E  01-0.10778855E-02  0.96830354E-06 
0.18713971E-09-0.22571089E-12  0.3(4128202  04  0.49370009E  00 
02  J  9/650  2.0  0.0  0.0  O.G  300.000  5000.000 

0.36219521E  01  0.73618256E-03-0.19652219E-06  0.36201556E-10-0.28945623E-14 
-0. 12019822E  04  0.36150942E  01  0.36255980E  01-0. 18782183E-02  0. 70554543E-05 
-0.67635071E-08  0.21555977E-11-0.10475225E  04  0.43052769E  01 
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FIGURE  5.  TYPICAL  THERMODYNAMIC  AND  REACTION  KINETICS  DATA  FILE  FOR  UMLAC 
AND  CREK  (LOGICAL  UNIT  15) 
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FIGURE  6.  (Continued) 
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